{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Multiple Kernel Learning"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "By Saurabh Mahindre - <a href=\"https://github.com/Saurabh7\">github.com/Saurabh7</a>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook is about multiple kernel learning in shogun. We will see how to construct a combined kernel, determine optimal kernel weights using MKL and use it for different types of [classification](http://en.wikipedia.org/wiki/Statistical_classification) and [novelty detection](http://en.wikipedia.org/wiki/Novelty_detection)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "1. [Introduction](#Introduction)\n",
      "2. [Mathematical formulation](#Mathematical-formulation-(skip if you just want code examples))\n",
      "3. [Using a Combined kernel](#Using-a-Combined-kernel)\n",
      "4. [Example: Toy Data](#Prediction-on-toy-data)\n",
      "  1. [Generating Kernel weights](#Generating-Kernel-weights)\n",
      "5. [Binary classification using MKL](#Binary-classification-using-MKL)\n",
      "6. [MKL for knowledge discovery](#MKL-for-knowledge-discovery)\n",
      "7. [Multiclass classification using MKL](#Multiclass-classification-using-MKL)\n",
      "8. [One-class classification using MKL](#One-class-classification-using-MKL)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "%matplotlib inline\n",
      "# import all shogun classes\n",
      "from modshogun import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Introduction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<em>Multiple kernel learning</em> (MKL) is about using a combined kernel i.e. a kernel consisting of a linear combination of arbitrary kernels over different domains. The coefficients or weights of the linear combination can be learned as well.\n",
      "\n",
      "[Kernel based methods](http://en.wikipedia.org/wiki/Kernel_methods) such as support vector machines (SVMs)  employ a so-called kernel function $k(x_{i},x_{j})$  which intuitively computes the similarity between two examples $x_{i}$ and $x_{j}$. </br>\n",
      "Selecting the kernel function\n",
      "$k()$  and it's parameters is an important issue in training. Kernels designed by humans usually capture one aspect of data. Choosing one kernel means to select exactly one such aspect. Which means combining such aspects is often better than selecting.\n",
      "\n",
      "In shogun the [CMKL](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKL.html) is the base class for MKL. We can do classifications:  [binary](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKLClassification.html), [one-class](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKLOneClass.html), [multiclass](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKLMulticlass.html) and regression too: [regression](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKLRegression.html)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Mathematical formulation (skip if you just want code examples)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "</br>In a SVM, defined as:\n",
      "$$f({\\bf x})=\\text{sign} \\left(\\sum_{i=0}^{N-1} \\alpha_i k({\\bf x}, {\\bf x_i})+b\\right)$$</br>\n",
      "where ${\\bf x_i},{i = 1,...,N}$ are labeled training examples ($y_i \\in {\u00b11}$).\n",
      "\n",
      "One could make a combination of kernels like:\n",
      "$${\\bf k}(x_i,x_j)=\\sum_{k=0}^{K} \\beta_k {\\bf k_k}(x_i, x_j)$$\n",
      "where $\\beta_k > 0$ and   $\\sum_{k=0}^{K} \\beta_k = 1$\n",
      "\n",
      "In the multiple kernel learning problem for binary classification one is given $N$ data points ($x_i, y_i$ )\n",
      "    ($y_i \\in {\u00b11}$), where $x_i$ is translated via $K$ mappings $\\phi_k(x) \\rightarrow R^{D_k} $, $k=1,...,K$ , from the input into $K$ feature spaces $(\\phi_1(x_i),...,\\phi_K(x_i))$ where $D_k$ denotes dimensionality of the $k$-th feature space.\n",
      "\n",
      "In MKL $\\alpha_i$,$\\beta$ and bias are determined by solving the following optimization program. For details see [1].\n",
      "\n",
      "$$\\mbox{min} \\hspace{4mm} \\gamma-\\sum_{i=1}^N\\alpha_i$$\n",
      "$$ \\mbox{w.r.t.} \\hspace{4mm}   \\gamma\\in R, \\alpha\\in R^N \\nonumber$$\n",
      "$$\\mbox {s.t.} \\hspace{4mm}  {\\bf 0}\\leq\\alpha\\leq{\\bf 1}C,\\;\\;\\sum_{i=1}^N \\alpha_i y_i=0 \\nonumber$$\n",
      "$$  {\\frac{1}{2}\\sum_{i,j=1}^N \\alpha_i \\alpha_j y_i y_j \\leq \\gamma},  \\forall k=1,\\ldots,K\\nonumber\\\\\n",
      "$$\n",
      "\n",
      "\n",
      "Here C is a pre-specified regularization parameter.\n",
      "Within shogun this optimization problem is solved using [semi-infinite programming](http://en.wikipedia.org/wiki/Semi-infinite_programming). For 1-norm MKL one of the two approaches described in [1] is used.\n",
      "The first approach (also called the wrapper algorithm) wraps around a single kernel SVMs, alternatingly solving for $\\alpha$ and $\\beta$. It is using a traditional SVM to generate new violated constraints and thus requires a single kernel SVM and any of the SVMs contained in shogun can be used. In the MKL step either a linear program is solved via [glpk](http://en.wikipedia.org/wiki/GNU_Linear_Programming_Kit) or cplex or analytically or a newton (for norms>1) step is performed.\n",
      "\n",
      "The second much faster but also more memory demanding approach performing interleaved optimization, is integrated into the chunking-based [SVMlight](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSVMLight.html).\n",
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Using a Combined kernel"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Shogun provides an easy way to make combination of kernels using the [CombinedKernel](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CCombinedKernel.html) class, to which we can append any [kernel](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKernel.html) from the many options shogun provides. It is especially useful to combine kernels working on different domains and to combine kernels looking at independent features and requires [CombinedFeatures](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CCombinedFeatures.html) to be used. Similarly the CombinedFeatures is used to combine a number of feature objects into a single CombinedFeatures object"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "kernel = CombinedKernel()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Prediction on toy data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In order to see the prediction capabilities, let us generate some data using the [GMM](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CGMM.html) class. The data is sampled by setting means ([GMM notebook](http://www.shogun-toolbox.org/static/notebook/current/GMM.html)) such that it sufficiently covers X-Y grid and is not too easy to classify."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "num=30;\n",
      "num_components=4\n",
      "means=zeros((num_components, 2))\n",
      "means[0]=[-1,1]\n",
      "means[1]=[2,-1.5]\n",
      "means[2]=[-1,-3]\n",
      "means[3]=[2,1]\n",
      "\n",
      "covs=array([[1.0,0.0],[0.0,1.0]])\n",
      "\n",
      "gmm=GMM(num_components)\n",
      "[gmm.set_nth_mean(means[i], i) for i in range(num_components)]\n",
      "[gmm.set_nth_cov(covs,i) for i in range(num_components)]\n",
      "gmm.set_coef(array([1.0,0.0,0.0,0.0]))\n",
      "xntr=array([gmm.sample() for i in xrange(num)]).T\n",
      "xnte=array([gmm.sample() for i in xrange(5000)]).T\n",
      "gmm.set_coef(array([0.0,1.0,0.0,0.0]))\n",
      "xntr1=array([gmm.sample() for i in xrange(num)]).T\n",
      "xnte1=array([gmm.sample() for i in xrange(5000)]).T\n",
      "gmm.set_coef(array([0.0,0.0,1.0,0.0]))\n",
      "xptr=array([gmm.sample() for i in xrange(num)]).T\n",
      "xpte=array([gmm.sample() for i in xrange(5000)]).T\n",
      "gmm.set_coef(array([0.0,0.0,0.0,1.0]))\n",
      "xptr1=array([gmm.sample() for i in xrange(num)]).T\n",
      "xpte1=array([gmm.sample() for i in xrange(5000)]).T\n",
      "traindata=concatenate((xntr,xntr1,xptr,xptr1), axis=1)\n",
      "trainlab=concatenate((-ones(2*num), ones(2*num)))\n",
      "\n",
      "testdata=concatenate((xnte,xnte1,xpte,xpte1), axis=1)\n",
      "testlab=concatenate((-ones(10000), ones(10000)))\n",
      "\n",
      "#convert to shogun features and generate labels for data\n",
      "feats_train=RealFeatures(traindata)  \n",
      "labels=BinaryLabels(trainlab)         "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "_=jet()\n",
      "figure(figsize(18,5))\n",
      "subplot(121)\n",
      "# plot train data\n",
      "_=scatter(traindata[0,:], traindata[1,:], c=trainlab, s=100)\n",
      "title('Toy data for classification')\n",
      "axis('equal')\n",
      "colors=[\"blue\",\"blue\",\"red\",\"red\"]\n",
      "# a tool for visualisation\n",
      "from matplotlib.patches import Ellipse\n",
      "def get_gaussian_ellipse_artist(mean, cov, nstd=1.96, color=\"red\", linewidth=3):\n",
      "    vals, vecs = eigh(cov)\n",
      "    order = vals.argsort()[::-1]\n",
      "    vals, vecs = vals[order], vecs[:, order]    \n",
      "    theta = numpy.degrees(arctan2(*vecs[:, 0][::-1]))\n",
      "    width, height = 2 * nstd * sqrt(vals)\n",
      "    e = Ellipse(xy=mean, width=width, height=height, angle=theta, \\\n",
      "               edgecolor=color, fill=False, linewidth=linewidth)\n",
      "    \n",
      "    return e\n",
      "for i in range(num_components):\n",
      "    gca().add_artist(get_gaussian_ellipse_artist(means[i], covs, color=colors[i]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Generating Kernel weights"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Just to help us visualize let's use two gaussian kernels ([CGaussianKernel](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CGaussianKernel.html)) with considerably different widths. As required in MKL, we need to append them to the Combined kernel. To generate the optimal weights (i.e $\\beta$s in the above equation), training of [MKL](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKLClassification.html) is required. This generates the weights as seen in this example."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "width0=0.5\n",
      "kernel0=GaussianKernel(feats_train, feats_train, width0)\n",
      "\n",
      "width1=25\n",
      "kernel1=GaussianKernel(feats_train, feats_train, width1)\n",
      "\n",
      "#combine kernels\n",
      "kernel.append_kernel(kernel0)   \n",
      "kernel.append_kernel(kernel1)\n",
      "kernel.init(feats_train, feats_train)\n",
      "\n",
      "mkl = MKLClassification()\n",
      "#set the norm, weights sum to 1.\n",
      "mkl.set_mkl_norm(1)  \n",
      "mkl.set_C(1, 1)\n",
      "mkl.set_kernel(kernel)\n",
      "mkl.set_labels(labels)\n",
      "\n",
      "#train to get weights\n",
      "mkl.train()    \n",
      "\n",
      "w=kernel.get_subkernel_weights()\n",
      "print w"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Binary classification using MKL"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now with the data ready and training done, we can do the binary classification. The weights generated can be intuitively understood. We will see that on plotting individual subkernels outputs and outputs of the MKL classification. To apply on test features, we need to reinitialize the kernel with `kernel.init` and pass the test features. After that it's just a matter of doing  `mkl.apply` to generate outputs. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size=100\n",
      "x1=linspace(-5, 5, size)\n",
      "x2=linspace(-5, 5, size)\n",
      "x, y=meshgrid(x1, x2)\n",
      "#Generate X-Y grid test data\n",
      "grid=RealFeatures(array((ravel(x), ravel(y))))\n",
      "\n",
      "kernel0t=GaussianKernel(feats_train, grid, width0)\n",
      "kernel1t=GaussianKernel(feats_train, grid, width1)\n",
      "\n",
      "kernelt=CombinedKernel()\n",
      "kernelt.append_kernel(kernel0t)\n",
      "kernelt.append_kernel(kernel1t)\n",
      "#initailize with test grid\n",
      "kernelt.init(feats_train, grid)\n",
      "\n",
      "mkl.set_kernel(kernelt)\n",
      "#prediction\n",
      "grid_out=mkl.apply()    \n",
      "\n",
      "z=grid_out.get_values().reshape((size, size))\n",
      "\n",
      "figure(figsize=(10,5))\n",
      "title(\"Classification using MKL\")\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To justify the weights, let's train and compare two subkernels with the MKL classification output. Training MKL classifier with a single kernel appended to a combined kernel makes no sense and is just like normal single kernel based classification, but let's do it for comparison."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "z=grid_out.get_labels().reshape((size, size))\n",
      "\n",
      "# MKL\n",
      "figure(figsize=(20,5))\n",
      "subplot(131, title=\"Multiple Kernels combined\")\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "\n",
      "comb_ker0=CombinedKernel()\n",
      "comb_ker0.append_kernel(kernel0)\n",
      "comb_ker0.init(feats_train, feats_train)\n",
      "mkl.set_kernel(comb_ker0)\n",
      "mkl.train()\n",
      "comb_ker0t=CombinedKernel()\n",
      "comb_ker0t.append_kernel(kernel0)\n",
      "comb_ker0t.init(feats_train, grid)\n",
      "mkl.set_kernel(comb_ker0t)\n",
      "out0=mkl.apply()\n",
      "\n",
      "# subkernel 1\n",
      "z=out0.get_labels().reshape((size, size)) \n",
      "subplot(132, title=\"Kernel 1\")\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "\n",
      "comb_ker1=CombinedKernel()\n",
      "comb_ker1.append_kernel(kernel1)\n",
      "comb_ker1.init(feats_train, feats_train)\n",
      "mkl.set_kernel(comb_ker1)\n",
      "mkl.train()\n",
      "comb_ker1t=CombinedKernel()\n",
      "comb_ker1t.append_kernel(kernel1)\n",
      "comb_ker1t.init(feats_train, grid)\n",
      "mkl.set_kernel(comb_ker1t)\n",
      "out1=mkl.apply()\n",
      "\n",
      "# subkernel 2\n",
      "z=out1.get_labels().reshape((size, size))  \n",
      "subplot(133, title=\"kernel 2\")\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As we can see the multiple kernel output seems just about right. Kernel 1 gives a sort of overfitting output while the kernel 2 seems not so accurate. The kernel weights are hence so adjusted to get a refined output. We can have a look at the errors by these subkernels to have more food for thought. Most of the time, the MKL error is lesser as it incorporates aspects of both kernels. One of them is strict while other is lenient, MKL finds a balance between those."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "kernelt.init(feats_train, RealFeatures(testdata))\n",
      "mkl.set_kernel(kernelt)\n",
      "out=mkl.apply()\n",
      "\n",
      "evaluator=ErrorRateMeasure()\n",
      "print \"Test error is %2.2f%% :MKL\" % (100*evaluator.evaluate(out,BinaryLabels(testlab)))\n",
      "\n",
      "\n",
      "comb_ker0t.init(feats_train,RealFeatures(testdata)) \n",
      "mkl.set_kernel(comb_ker0t)\n",
      "out=mkl.apply()\n",
      "\n",
      "evaluator=ErrorRateMeasure()\n",
      "print \"Test error is %2.2f%% :Subkernel1\"% (100*evaluator.evaluate(out,BinaryLabels(testlab)))\n",
      "\n",
      "comb_ker1t.init(feats_train, RealFeatures(testdata))\n",
      "mkl.set_kernel(comb_ker1t)\n",
      "out=mkl.apply()\n",
      "\n",
      "evaluator=ErrorRateMeasure()\n",
      "print \"Test error is %2.2f%% :subkernel2\" % (100*evaluator.evaluate(out,BinaryLabels(testlab)))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "MKL for knowledge discovery"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "MKL can recover information about the problem at hand. Let us see this with a binary classification problem. The task is to separate two concentric classes shaped like circles. By varying the distance between the boundary of the circles we can control the separability of the problem. Starting with an almost non-separable scenario, the data quickly becomes separable as the distance between the circles increases."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def circle(x, radius, neg):\n",
      "        y=sqrt(square(radius)-square(x))\n",
      "        if neg:\n",
      "            return[x, -y]\n",
      "        else:\n",
      "            return [x,y]\n",
      "        \n",
      "def get_circle(radius):\n",
      "    neg=False\n",
      "    range0=linspace(-radius,radius,100)\n",
      "    pos_a=array([circle(i, radius, neg) for i in range0]).T\n",
      "    neg=True\n",
      "    neg_a=array([circle(i, radius, neg) for i in range0]).T\n",
      "    c=concatenate((neg_a,pos_a), axis=1)\n",
      "    return c\n",
      "\n",
      "def get_data(r1, r2):\n",
      "    c1=get_circle(r1)\n",
      "    c2=get_circle(r2)\n",
      "    c=concatenate((c1, c2), axis=1)\n",
      "    feats_tr=RealFeatures(c)\n",
      "    return c, feats_tr\n",
      "\n",
      "l=concatenate((-ones(200),ones(200)))\n",
      "lab=BinaryLabels(l)\n",
      "\n",
      "#get two circles with radius 2 and 4\n",
      "c, feats_tr=get_data(2,4)\n",
      "c1, feats_tr1=get_data(2,3)\n",
      "_=gray()\n",
      "figure(figsize=(10,5))\n",
      "subplot(121)\n",
      "title(\"Circles with different separation\")\n",
      "p=scatter(c[0,:], c[1,:], c=lab)\n",
      "subplot(122)\n",
      "q=scatter(c1[0,:], c1[1,:], c=lab)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These are the type of circles we want to distinguish between. We can try classification with a constant separation between the circles first."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def train_mkl(circles, feats_tr):\n",
      "    #Four kernels with different widths \n",
      "    kernel0=GaussianKernel(feats_tr, feats_tr, 1)  \n",
      "    kernel1=GaussianKernel(feats_tr, feats_tr, 5)\n",
      "    kernel2=GaussianKernel(feats_tr, feats_tr, 7)\n",
      "    kernel3=GaussianKernel(feats_tr, feats_tr, 10)\n",
      "    kernel = CombinedKernel()\n",
      "    kernel.append_kernel(kernel0)\n",
      "    kernel.append_kernel(kernel1)\n",
      "    kernel.append_kernel(kernel2)\n",
      "    kernel.append_kernel(kernel3)\n",
      "    \n",
      "    kernel.init(feats_tr, feats_tr)\n",
      "    mkl = MKLClassification()\n",
      "    mkl.set_mkl_norm(1)\n",
      "    mkl.set_C(1, 1)\n",
      "    mkl.set_kernel(kernel)\n",
      "    mkl.set_labels(lab)\n",
      "    \n",
      "    mkl.train()\n",
      "    \n",
      "    w=kernel.get_subkernel_weights()\n",
      "    return w, mkl\n",
      "\n",
      "def test_mkl(mkl, grid):\n",
      "    kernel0t=GaussianKernel(feats_tr, grid, 1)\n",
      "    kernel1t=GaussianKernel(feats_tr, grid, 5)\n",
      "    kernel2t=GaussianKernel(feats_tr, grid, 7)\n",
      "    kernel3t=GaussianKernel(feats_tr, grid, 10)\n",
      "    kernelt = CombinedKernel()\n",
      "    kernelt.append_kernel(kernel0t)\n",
      "    kernelt.append_kernel(kernel1t)\n",
      "    kernelt.append_kernel(kernel2t)\n",
      "    kernelt.append_kernel(kernel3t)\n",
      "    kernelt.init(feats_tr, grid)\n",
      "    mkl.set_kernel(kernelt)\n",
      "    out=mkl.apply()\n",
      "    return out\n",
      "\n",
      "size=50\n",
      "x1=linspace(-10, 10, size)\n",
      "x2=linspace(-10, 10, size)\n",
      "x, y=meshgrid(x1, x2)\n",
      "grid=RealFeatures(array((ravel(x), ravel(y))))\n",
      "\n",
      "\n",
      "w, mkl=train_mkl(c, feats_tr)\n",
      "print w\n",
      "out=test_mkl(mkl,grid)\n",
      "\n",
      "z=out.get_values().reshape((size, size))\n",
      "\n",
      "figure(figsize=(5,5))\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "title('classification with constant separation')\n",
      "_=colorbar(c)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As we can see the MKL classifier classifies them as expected. Now let's vary the separation and see how it affects the weights.The choice of the kernel width of the Gaussian kernel used for classification is expected to depend on the separation distance of the learning problem. An increased distance between the circles will correspond to a larger optimal kernel width. This effect should be visible in the results of the MKL, where we used MKL-SVMs with four kernels with different widths (1,5,7,10). "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "range1=linspace(5.5,7.5,50)\n",
      "x=linspace(1.5,3.5,50)\n",
      "temp=[]\n",
      "\n",
      "for i in range1:\n",
      "    #vary separation between circles\n",
      "    c, feats=get_data(4,i) \n",
      "    w, mkl=train_mkl(c, feats)\n",
      "    temp.append(w)\n",
      "y=array([temp[i] for i in range(0,50)]).T\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure(figsize=(20,5))\n",
      "_=plot(x, y[0,:], color='k', linewidth=2)\n",
      "_=plot(x, y[1,:], color='r', linewidth=2)\n",
      "_=plot(x, y[2,:], color='g', linewidth=2)\n",
      "_=plot(x, y[3,:], color='y', linewidth=2)\n",
      "title(\"Comparison between kernel widths and weights\")\n",
      "ylabel(\"Weight\")\n",
      "xlabel(\"Distance between circles\")\n",
      "_=legend([\"1\",\"5\",\"7\",\"10\"])\n",
      " "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the above plot we see the kernel weightings obtained for the four kernels. Every line shows one weighting. The courses of the kernel weightings reflect the development of the learning problem: as long as the problem is difficult the best separation can be obtained when using the kernel with smallest width. The low width kernel looses importance when the distance between the circle increases and larger kernel widths obtain a larger weight in MKL. Increasing the distance between the circles, kernels with greater widths are used. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Multiclass classification using MKL"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "MKL can be used for multiclass classification using the [MKLMulticlass](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKLMulticlass.html) class. It is based on the GMNPSVM Multiclass SVM. Its termination criterion is set by  `set_mkl_epsilon(float64_t eps )` and the maximal number of MKL iterations is set by `set_max_num_mkliters(int32_t maxnum)`. The epsilon termination criterion is the L2 norm between the current MKL weights and their counterpart from the previous iteration. We set it to 0.001 as we want pretty accurate weights.\n",
      "\n",
      "To see this in action let us compare it to the normal [GMNPSVM](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CGMNPSVM.html) example as in the [KNN notebook](http://www.shogun-toolbox.org/static/notebook/current/KNN.html#Comparison-to-Multiclass-Support-Vector-Machines), just to see how MKL fares in object recognition. We use the  [USPS digit recognition dataset](http://www.gaussianprocess.org/gpml/data/)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.io import loadmat, savemat\n",
      "from os       import path, sep\n",
      "\n",
      "mat  = loadmat(sep.join(['..','..','..','data','multiclass', 'usps.mat']))\n",
      "Xall = mat['data']\n",
      "Yall = array(mat['label'].squeeze(), dtype=double)\n",
      "\n",
      "# map from 1..10 to 0..9, since shogun\n",
      "# requires multiclass labels to be\n",
      "# 0, 1, ..., K-1\n",
      "Yall = Yall - 1\n",
      "\n",
      "random.seed(0)\n",
      "\n",
      "subset = random.permutation(len(Yall))\n",
      "\n",
      "#get first 1000 examples\n",
      "Xtrain = Xall[:, subset[:1000]]\n",
      "Ytrain = Yall[subset[:1000]]\n",
      "\n",
      "Nsplit = 2\n",
      "all_ks = range(1, 21)\n",
      "\n",
      "print Xall.shape\n",
      "print Xtrain.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's plot five of the  examples to get a feel of the dataset."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def plot_example(dat, lab):\n",
      "    for i in xrange(5):\n",
      "        ax=subplot(1,5,i+1)\n",
      "        title(int(lab[i]))\n",
      "        ax.imshow(dat[:,i].reshape((16,16)), interpolation='nearest')\n",
      "        ax.set_xticks([])\n",
      "        ax.set_yticks([])\n",
      "        \n",
      "        \n",
      "_=figure(figsize=(17,6))\n",
      "gray()\n",
      "plot_example(Xtrain, Ytrain)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We combine a [Gaussian kernel](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CGaussianKernel.html) and a [PolyKernel](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CPolyKernel.html). To test, examples not included in training data are used.\n",
      "\n",
      "This is just a demonstration but we can see here how MKL is working behind the scene. What we have is two kernels with significantly different properties. The gaussian kernel defines a function space that is a lot larger than that of the linear kernel or the polynomial kernel. The gaussian kernel has a low width, so it will be able to represent more and more complex relationships between the training data. But it requires enough data to train on. The number of training examples here is 1000, which seems a bit less as total examples are 10000. We hope the polynomial kernel can counter this problem, since it will fit the polynomial for you using a lot less data than the squared exponential. The kernel weights are printed below to add some insight."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# MKL training and output\n",
      "labels = MulticlassLabels(Ytrain)\n",
      "feats  = RealFeatures(Xtrain)\n",
      "#get test data from 5500 onwards\n",
      "Xrem=Xall[:,subset[5500:]]\n",
      "Yrem=Yall[subset[5500:]]\n",
      "\n",
      "#test features not used in training\n",
      "feats_rem=RealFeatures(Xrem)             \n",
      "labels_rem=MulticlassLabels(Yrem)\n",
      "\n",
      "kernel = CombinedKernel()\n",
      "feats_train = CombinedFeatures()\n",
      "feats_test = CombinedFeatures()\n",
      "\n",
      "#append gaussian kernel\n",
      "subkernel = GaussianKernel(10,15)        \n",
      "feats_train.append_feature_obj(feats)\n",
      "feats_test.append_feature_obj(feats_rem)\n",
      "kernel.append_kernel(subkernel)\n",
      "\n",
      "#append PolyKernel\n",
      "feats  = RealFeatures(Xtrain)\n",
      "subkernel = PolyKernel(10,2)            \n",
      "feats_train.append_feature_obj(feats)\n",
      "feats_test.append_feature_obj(feats_rem)\n",
      "kernel.append_kernel(subkernel)\n",
      "\n",
      "kernel.init(feats_train, feats_train)\n",
      "\n",
      "mkl = MKLMulticlass(1.2, kernel, labels)\n",
      "\n",
      "mkl.set_epsilon(1e-2)\n",
      "mkl.set_mkl_epsilon(0.001)\n",
      "mkl.set_mkl_norm(1)\n",
      "\n",
      "mkl.train()\n",
      "\n",
      "#initialize with test features\n",
      "kernel.init(feats_train, feats_test)     \n",
      "\n",
      "out =  mkl.apply()\n",
      "evaluator = MulticlassAccuracy()\n",
      "accuracy = evaluator.evaluate(out, labels_rem)\n",
      "print \"Accuracy = %2.2f%%\" % (100*accuracy)\n",
      "\n",
      "idx=where(out.get_labels() != Yrem)[0]\n",
      "Xbad=Xrem[:,idx]\n",
      "Ybad=Yrem[idx]\n",
      "_=figure(figsize=(17,6))\n",
      "gray()\n",
      "plot_example(Xbad, Ybad)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w=kernel.get_subkernel_weights()\n",
      "print w"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Single kernel:PolyKernel\n",
      "C=1\n",
      "\n",
      "pk=PolyKernel(10,2) \n",
      "\n",
      "svm=GMNPSVM(C, pk, labels)\n",
      "_=svm.train(feats)\n",
      "out=svm.apply(feats_rem)\n",
      "evaluator = MulticlassAccuracy()\n",
      "accuracy = evaluator.evaluate(out, labels_rem)\n",
      "\n",
      "print \"Accuracy = %2.2f%%\" % (100*accuracy)\n",
      "\n",
      "idx=np.where(out.get_labels() != Yrem)[0]\n",
      "Xbad=Xrem[:,idx]\n",
      "Ybad=Yrem[idx]\n",
      "_=figure(figsize=(17,6))\n",
      "gray()\n",
      "plot_example(Xbad, Ybad)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Single Kernel:Gaussian kernel\n",
      "width=15\n",
      "C=1\n",
      "\n",
      "gk=GaussianKernel()\n",
      "gk.set_width(width)\n",
      "\n",
      "svm=GMNPSVM(C, gk, labels)\n",
      "_=svm.train(feats)\n",
      "out=svm.apply(feats_rem)\n",
      "evaluator = MulticlassAccuracy()\n",
      "accuracy = evaluator.evaluate(out, labels_rem)\n",
      "\n",
      "print \"Accuracy = %2.2f%%\" % (100*accuracy)\n",
      "\n",
      "idx=np.where(out.get_labels() != Yrem)[0]\n",
      "Xbad=Xrem[:,idx]\n",
      "Ybad=Yrem[idx]\n",
      "_=figure(figsize=(17,6))\n",
      "gray()\n",
      "plot_example(Xbad, Ybad)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The misclassified examples are surely pretty tough to predict. As seen from the accuracy MKL seems to work a shade better in the case. One could try this out with more and different types of kernels too."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "One-class classification using MKL"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[One-class classification](http://en.wikipedia.org/wiki/One-class_classification) can be done using MKL in shogun. This is demonstrated in the following simple example using [CMKLOneClass](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMKLOneClass.html). We will see how abnormal data is detected. This is also known as novelty detection. Below we generate some toy data and initialize combined kernels and features."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = -0.3 * random.randn(100,2)\n",
      "traindata=r_[X + 2, X - 2].T\n",
      "\n",
      "X = -0.3 * random.randn(20, 2)\n",
      "testdata = r_[X + 2, X - 2].T\n",
      "\n",
      "trainlab=concatenate((ones(99),-ones(1)))\n",
      "#convert to shogun features and generate labels for data\n",
      "feats=RealFeatures(traindata) \n",
      "labels=BinaryLabels(trainlab)         "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xx, yy = meshgrid(linspace(-5, 5, 500), linspace(-5, 5, 500))\n",
      "grid=RealFeatures(array((ravel(xx), ravel(yy))))\n",
      "#test features\n",
      "feats_t=RealFeatures(testdata)   \n",
      "x_out=(random.uniform(low=-4, high=4, size=(20, 2))).T\n",
      "feats_out=RealFeatures(x_out)\n",
      "\n",
      "kernel=CombinedKernel()\n",
      "feats_train=CombinedFeatures()\n",
      "feats_test=CombinedFeatures()\n",
      "feats_test_out=CombinedFeatures()\n",
      "feats_grid=CombinedFeatures()\n",
      "\n",
      "#append gaussian kernel\n",
      "subkernel=GaussianKernel(10,8)        \n",
      "feats_train.append_feature_obj(feats)\n",
      "feats_test.append_feature_obj(feats_t)\n",
      "feats_test_out.append_feature_obj(feats_out)\n",
      "feats_grid.append_feature_obj(grid)\n",
      "kernel.append_kernel(subkernel)\n",
      "\n",
      "#append PolyKernel\n",
      "feats  = RealFeatures(traindata)\n",
      "subkernel = PolyKernel(10,3)            \n",
      "feats_train.append_feature_obj(feats)\n",
      "feats_test.append_feature_obj(feats_t)\n",
      "feats_test_out.append_feature_obj(feats_out)\n",
      "feats_grid.append_feature_obj(grid)\n",
      "kernel.append_kernel(subkernel)\n",
      "\n",
      "kernel.init(feats_train, feats_train)\n",
      "\n",
      "mkl = MKLOneClass()\n",
      "mkl.set_kernel(kernel)\n",
      "mkl.set_labels(labels)\n",
      "mkl.set_interleaved_optimization_enabled(False)\n",
      "\n",
      "mkl.set_epsilon(1e-2)\n",
      "mkl.set_mkl_epsilon(0.1)\n",
      "mkl.set_mkl_norm(1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that everything is initialized, let's see MKLOneclass in action by applying it on the test data and on the X-Y grid."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mkl.train()\n",
      "print \"Weights:\"\n",
      "w=kernel.get_subkernel_weights()\n",
      "print w\n",
      "\n",
      "#initialize with test features\n",
      "kernel.init(feats_train, feats_test)     \n",
      "normal_out =  mkl.apply()\n",
      "\n",
      "#test on abnormally generated data\n",
      "kernel.init(feats_train, feats_test_out)     \n",
      "abnormal_out =  mkl.apply()\n",
      "\n",
      "#test on X-Y grid\n",
      "kernel.init(feats_train, feats_grid)\n",
      "grid_out=mkl.apply()\n",
      "z=grid_out.get_values().reshape((500,500))\n",
      "z_lab=grid_out.get_labels().reshape((500,500))\n",
      "\n",
      "a=abnormal_out.get_labels()\n",
      "n=normal_out.get_labels()\n",
      "\n",
      "#check for normal and abnormal classified data\n",
      "idx=where(normal_out.get_labels() != 1)[0]\n",
      "abnormal=testdata[:,idx]\n",
      "\n",
      "idx=where(normal_out.get_labels() == 1)[0]\n",
      "normal=testdata[:,idx]\n",
      "\n",
      "figure(figsize(15,6))\n",
      "pl =subplot(121)\n",
      "title(\"One-class classification using MKL\")\n",
      "_=pink()\n",
      "c=pcolor(xx, yy, z)\n",
      "_=contour(xx, yy, z_lab, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "p1=pl.scatter(traindata[0, :], traindata[1,:], c=labels,cmap=gray(), s=100)\n",
      "p2=pl.scatter(normal[0,:], normal[1,:], c=\"red\", s=100)\n",
      "p3=pl.scatter(abnormal[0,:], abnormal[1,:], c=\"blue\", s=100)\n",
      "p4=pl.scatter(x_out[0,:], x_out[1,:], c=a, cmap=jet(),  s=100)\n",
      "_=pl.legend((p1, p2, p3), [\"Training samples\", \"normal samples\", \"abnormal samples\"], loc=2)\n",
      "\n",
      "\n",
      "\n",
      "subplot(122)\n",
      "c=pcolor(xx, yy, z)\n",
      "title(\"One-class classification output\")\n",
      "_=gray()\n",
      "_=contour(xx, yy, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "MKL one-class classification will give you a bit more flexibility compared to normal classifications. The kernel weights are expected to be more or less similar here since the training data is not overly complicated or too easy, which means both the gaussian and polynomial kernel will be involved. If you don't know the nature of the training data and lot of features are invoved, you could easily use kernels with much different properties and benefit from their combination."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "References:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[1] Soeren Sonnenburg, Gunnar Raetsch, Christin Schaefer, and Bernhard Schoelkopf. Large Scale Multiple Kernel Learning. Journal of Machine Learning Research, 7:1531-1565, July 2006.\n",
      "\n",
      "[2]F. R. Bach, G. R. G. Lanckriet, and M. I. Jordan. Multiple kernel learning, conic duality, and\n",
      "the SMO algorithm. In C. E. Brodley, editor, Twenty-first international conference on Machine\n",
      "learning. ACM, 2004\n",
      "\n",
      "[3] Kernel Methods for Object Recognition , Christoph H. Lampert"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}
